Modelando Preços de Venda de Imóveis com Modelo Linear Normal
Elton Dantas de Oliveira Mesquita
Pacotes
library(tidyverse)
library(plotly)
library(GGally)
library(gridExtra)
library(ggfortify)
library(hnp)
library(nortest)Dados
Exemplo disponível em: Paula, G. A. (2013). Modelos de Regressão com Apoio Computacional. São Paulo, SP: IME-USP. (Exercício 23, página 112). Dados disponíveis em: https://www.ime.usp.br/~giapaula/textoregressao.htm
Neste trabalho, vamos modelar o preço de venda de imóveis a partir de dados relativos a uma amostra de 27 imóveis. As variáveis desse conjunto de dados são:
imposto: imposto do imóvel (em US$ 100);
areaT: área do terreno (em 1.000 pés quadrados);
areaC: área construída (em 1.000 pés quadrados);
idade: idade da residência (em anos);
preco: preço de venda do imóvel (em US$ 1.000).
Sendo assim, o objetivo dessa modelo é explicar a variável preço de venda a partir das variáveis área do terreno, área construída e idade.
dados = read.table("imoveis.dat")
colnames(dados) = c("imposto", "areaT", "areaC", "idade", "preco")
attach(dados)
dados## imposto areaT areaC idade preco
## 1 4.9176 3.4720 0.998 42 25.9
## 2 5.0208 3.5310 1.500 62 29.5
## 3 4.5429 2.2750 1.175 40 27.9
## 4 4.5573 4.0500 1.232 54 25.9
## 5 5.0597 4.4550 1.121 42 29.9
## 6 3.8910 4.4550 0.988 56 29.9
## 7 5.8980 5.8500 1.240 51 30.9
## 8 5.6039 9.5200 1.501 32 28.9
## 9 15.4202 9.8000 3.420 42 84.9
## 10 14.4598 12.8000 3.000 14 82.9
## 11 5.8282 6.4350 1.225 32 35.9
## 12 5.3003 4.9883 1.552 30 31.5
## 13 6.2712 5.5200 0.975 30 31.0
## 14 5.9592 6.6660 1.121 32 30.9
## 15 5.0500 5.0000 1.020 46 30.0
## 16 8.2464 5.1500 1.664 50 36.9
## 17 6.6969 6.9020 1.488 22 41.9
## 18 7.7841 7.1020 1.376 17 40.5
## 19 9.0384 7.8000 1.500 23 43.9
## 20 5.9894 5.5200 1.256 40 37.5
## 21 7.5422 4.0000 1.690 22 37.9
## 22 8.7951 9.8900 1.820 50 44.5
## 23 6.0931 6.7265 1.652 44 37.9
## 24 8.3607 9.1500 1.777 48 38.9
## 25 8.1400 8.0000 1.504 3 36.9
## 26 9.1416 7.3262 1.831 31 45.8
## 27 12.0000 5.0000 1.200 30 41.0Análise Descritiva
Box-plots
O gráfico interativo a baixo nos mostra algumas medidas resumo para cada uma das variáveis.
plot_ly(type = "box") %>%
add_trace(y = imposto, name = "imposto") %>%
add_trace(y = areaT, name = "área do terreno") %>%
add_trace(y = areaC, name = "área construída") %>%
add_trace(y = idade, name = "idade do imóvel") %>%
add_trace(y = preco, name = "preço de venda")Densidades, Dispersões e Correlações
g = ggpairs(dados, aes(color = I("slategray"), fill = I("slategray")),
lower = list(continuous = wrap("smooth",col="black")),
diag=list(continuous = wrap("densityDiag",alpha=0.5,size=1)))
ggplotly(g)Nos gráficos descritivos acima, percebemos que todas as variáveis possuem distribuições assimétricas e que a variável preço possui altas correlações positivas com as demais, exceto com a variável idade, em que a correlação é fraca e negativa. Em outras palavras, quanto maiores os valores de imposto, área do terreno e área construída, maior o valor do preço de venda. Também identificamos a presença de outliers, mais evidentes nos gráficos de box-plot.
Ajuste 1
Nesse primeiro ajuste, consideraremos para o nosso modelo todas as variáveis disponíveis e o chamaremos de modelo completo. Partiremos do modelo normal linear.
\[ preço = \beta_0 + \beta_1*imposto + \beta_2*areaC + \beta_3*areaT + \beta_4*idade + \epsilon \]
Modelo Completo
model1 = lm(preco~.,dados)
summary(model1)##
## Call:
## lm(formula = preco ~ ., data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9685 -2.7152 0.2663 2.1374 6.9872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.43587 4.09233 0.595 0.55776
## imposto 2.07823 0.55296 3.758 0.00109 **
## areaT 0.23238 0.50658 0.459 0.65094
## areaC 13.97381 2.90663 4.808 8.4e-05 ***
## idade -0.04376 0.06628 -0.660 0.51594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.077 on 22 degrees of freedom
## Multiple R-squared: 0.9313, Adjusted R-squared: 0.9188
## F-statistic: 74.56 on 4 and 22 DF, p-value: 1.809e-12Verificamos que, para o modelo completo, as variáveis imposto e área construída são significativas e obtivemos um \(R^2\) muito alto, indicando um possível bom ajuste aos dados amostrais. Porém, o \(R^2\) por si só não nos garante que o modelo seja o mais adequado.
A seguir, faremos a análise de diagnóstico para verificar a adequação do modelo.
Diagnóstico 1
Normalidade e Envelope
Traçaremos agora os gráficos de histograma e envelope, para avaliar o ajuste do modelo normal.
hist(model1$residuals, probability = TRUE, main = 'Resíduos do Ajuste',
xlab = 'Resíduos', ylab = 'Frequência', col = "gray",)
lines(density(model1$residuals), col = "slategray", lwd = 4)
hnp(model1, main = 'Gráfico Normal com Envelope',
xlab = 'Quantis Teóricos',
ylab = 'Resíduos',
pch = 16)## Gaussian model (lm object)
Testes de Normalidade dos Resíduos
Para os testes a seguir, consideramos
\(H_0\): Os resíduos têm distribuição normal.
# Kolmogorov - Smirnov
ks.test(model1$residuals,"pnorm")##
## One-sample Kolmogorov-Smirnov test
##
## data: model1$residuals
## D = 0.34874, p-value = 0.001937
## alternative hypothesis: two-sided# Lilliefors
lillie.test(model1$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: model1$residuals
## D = 0.069288, p-value = 0.9841# Cramer-von-Mises
cvm.test(model1$residuals)##
## Cramer-von Mises normality test
##
## data: model1$residuals
## W = 0.021868, p-value = 0.9451# Shapiro-Wilk
shapiro.test(model1$residuals)##
## Shapiro-Wilk normality test
##
## data: model1$residuals
## W = 0.98111, p-value = 0.8865# Shepiro-Francia
sf.test(model1$residuals)##
## Shapiro-Francia normality test
##
## data: model1$residuals
## W = 0.98974, p-value = 0.9768# Anderson-Darling
ad.test(model1$residuals)##
## Anderson-Darling normality test
##
## data: model1$residuals
## A = 0.15285, p-value = 0.9528Ao nível de 5% de significância, não rejeitamos a hipótese de que os resíduos tenham distribuição normal.
Análise de Resíduos
plot(lm(model1, dados), which = 1:4, pch = 16,lwd = 3,
caption = c("Resíduos x Valores ajustados",
"Normal Q-Q",
"Resíduos padronizados x Valores ajustados",
"Distância de Cook",
"Resíduos x Alavancagem"))
No gráfico de Resíduos x Valores Ajustados, verifica-se uma mudança na faixa de variação dos pontos, isto é, indícios de heterocedasticidade. Há também outliers
Medidas de Influência
DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo (apontam os pontos influentes);
DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores preditos/ajustados do modelo (apontam os pontos aberrantes);
COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo (aponta os pontos alavanca);
HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados.
inf = influence.measures(model1)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ ., data = dados) :
##
## dfb.1_ dfb.imps dfb.areT dfb.areC dfb.idad dffit cov.r cook.d hat
## 8 -0.01 0.73 -1.01_* -0.05 0.08 -1.29 0.63 0.28 0.26
## 9 -0.49 0.15 -0.28 0.47 0.18 0.97 2.58_* 0.19 0.58_*
## 10 -0.44 -0.02 0.40 0.44 -0.38 1.62_* 0.82 0.46 0.38
## 27 0.11 -1.91_* 0.47 1.41_* -0.36 -2.02_* 2.71_* 0.79 0.69_*# A soma da diagonal da matriz H é o número de parâmetros do modelo
# sum(inf$infmat[,"hat"])A partir da tabela acima, temos que:
Observação 8: ponto influente;
Observação 9: ponto de alavanca;
Observação 10: ponto aberrante;
Observação 27: ponto influente, aberrante e de alavanca.
Ajuste 2
A seguir, tentaremos melhorar o ajuste anterior, selecionando variáveis através da técnica AIC. Mas, além de todas as variáveis do modelo completo, consideraremos também as suas interações.
Seleção de Variáveis (AIC)
model = lm(preco~(.)^2,dados)
MASS::stepAIC(model)## Start: AIC=70.09
## preco ~ (imposto + areaT + areaC + idade)^2
##
## Df Sum of Sq RSS AIC
## - imposto:idade 1 0.232 160.53 68.132
## - areaC:idade 1 6.109 166.41 69.103
## - areaT:idade 1 6.565 166.87 69.177
## <none> 160.30 70.093
## - areaT:areaC 1 14.674 174.98 70.458
## - imposto:areaT 1 14.764 175.07 70.472
## - imposto:areaC 1 53.086 213.39 75.817
##
## Step: AIC=68.13
## preco ~ imposto + areaT + areaC + idade + imposto:areaT + imposto:areaC +
## areaT:areaC + areaT:idade + areaC:idade
##
## Df Sum of Sq RSS AIC
## - areaT:idade 1 10.974 171.51 67.918
## <none> 160.53 68.132
## - areaT:areaC 1 16.309 176.84 68.745
## - imposto:areaT 1 17.096 177.63 68.865
## - areaC:idade 1 18.031 178.57 69.006
## - imposto:areaC 1 53.024 213.56 73.838
##
## Step: AIC=67.92
## preco ~ imposto + areaT + areaC + idade + imposto:areaT + imposto:areaC +
## areaT:areaC + areaC:idade
##
## Df Sum of Sq RSS AIC
## - areaC:idade 1 8.891 180.40 67.282
## - imposto:areaT 1 10.456 181.96 67.515
## - areaT:areaC 1 10.934 182.44 67.586
## <none> 171.51 67.918
## - imposto:areaC 1 48.083 219.59 72.590
##
## Step: AIC=67.28
## preco ~ imposto + areaT + areaC + idade + imposto:areaT + imposto:areaC +
## areaT:areaC
##
## Df Sum of Sq RSS AIC
## - areaT:areaC 1 4.907 185.31 66.007
## - idade 1 10.805 191.20 66.853
## - imposto:areaT 1 13.674 194.07 67.255
## <none> 180.40 67.282
## - imposto:areaC 1 42.514 222.91 70.996
##
## Step: AIC=66.01
## preco ~ imposto + areaT + areaC + idade + imposto:areaT + imposto:areaC
##
## Df Sum of Sq RSS AIC
## - idade 1 14.233 199.54 66.005
## <none> 185.31 66.007
## - imposto:areaT 1 17.866 203.17 66.492
## - imposto:areaC 1 40.523 225.83 69.346
##
## Step: AIC=66
## preco ~ imposto + areaT + areaC + imposto:areaT + imposto:areaC
##
## Df Sum of Sq RSS AIC
## <none> 199.54 66.005
## - imposto:areaT 1 21.923 221.46 66.819
## - imposto:areaC 1 33.793 233.33 68.229##
## Call:
## lm(formula = preco ~ imposto + areaT + areaC + imposto:areaT +
## imposto:areaC, data = dados)
##
## Coefficients:
## (Intercept) imposto areaT areaC imposto:areaT
## 25.4793 -0.1035 -0.7113 -0.5733 0.1798
## imposto:areaC
## 0.8215Inicialmente, verificamos que imposto e área construída eram as únicas variáveis significativas para o modelo e elas também foram selecionadas pela técnica AIC. Considerando as interações das variáveis do modelo completo, a seleção por AIC nos sugere também a variável imposto:areaT, que relaciona o imposto à área do terreno. Chamaremos esse ajuste de modelo reduzido e repetiremos as técnicas de diagnóstico realizadas no modelo completo.
Modelo Reduzido
model2 = lm(preco~imposto+areaT+areaC+imposto:areaT,dados)
summary(model2)##
## Call:
## lm(formula = preco ~ imposto + areaT + areaC + imposto:areaT,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0536 -2.2749 -0.0292 2.0684 6.2473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.63620 5.94598 3.471 0.00217 **
## imposto 0.24245 0.68136 0.356 0.72535
## areaT -1.95662 0.73976 -2.645 0.01479 *
## areaC 7.76323 2.66836 2.909 0.00813 **
## imposto:areaT 0.33200 0.09153 3.627 0.00149 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.257 on 22 degrees of freedom
## Multiple R-squared: 0.9562, Adjusted R-squared: 0.9482
## F-statistic: 120 on 4 and 22 DF, p-value: 1.324e-14Diagnóstico 2
Normalidade e Envelope
hist(model2$residuals, probability = TRUE, main = 'Resíduos do Ajuste',
xlab = 'Resíduos', ylab = 'Frequência', col = "gray",)
lines(density(model$residuals), col = "slategray", lwd = 4)
hnp(model2, main = 'Gráfico Normal com Envelope',
xlab = 'Quantis Teóricos',
ylab = 'Resíduos',
pch = 16)## Gaussian model (lm object)
Testes de Normalidade dos Resíduos
# Kolmogorov - Smirnov
ks.test(model2$residuals,"pnorm")##
## One-sample Kolmogorov-Smirnov test
##
## data: model2$residuals
## D = 0.34327, p-value = 0.002409
## alternative hypothesis: two-sided# Lilliefors
lillie.test(model2$residuals)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: model2$residuals
## D = 0.11015, p-value = 0.5465# Cramer-von-Mises
cvm.test(model2$residuals)##
## Cramer-von Mises normality test
##
## data: model2$residuals
## W = 0.041833, p-value = 0.6335# Shapiro-Wilk
shapiro.test(model2$residuals)##
## Shapiro-Wilk normality test
##
## data: model2$residuals
## W = 0.96969, p-value = 0.5937# Shepiro-Francia
sf.test(model2$residuals)##
## Shapiro-Francia normality test
##
## data: model2$residuals
## W = 0.97702, p-value = 0.7# Anderson-Darling
ad.test(model2$residuals)##
## Anderson-Darling normality test
##
## data: model2$residuals
## A = 0.27501, p-value = 0.6331Análise de Resíduos
plot(lm(model2, dados), which = 1:4, pch = 16,lwd = 3,
caption = c("Resíduos x Valores ajustados",
"Normal Q-Q",
"Resíduos padronizados x Valores ajustados",
"Distância de Cook",
"Resíduos x Alavancagem"))
Medidas de Influência
DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo (apontam os pontos influentes);
DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores preditos/ajustados do modelo (apontam os pontos aberrantes);
COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo (aponta os pontos alavanca);
HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados.
inf = influence.measures(model2)
summary(inf)## Potentially influential observations of
## lm(formula = preco ~ imposto + areaT + areaC + imposto:areaT, data = dados) :
##
## dfb.1_ dfb.imps dfb.areT dfb.areC dfb.im:T dffit cov.r cook.d hat
## 9 -0.53 0.26 -0.13 0.92 -0.22 1.65_* 1.85_* 0.52 0.57_*
## 10 -0.52 0.56 0.51 0.31 -0.75 -1.03 4.88_* 0.22 0.76_*
## 27 0.14 -1.11_* 0.06 0.77 0.23 -1.61_* 3.02_* 0.52 0.68_*# A soma da diagonal da matriz H é o número de parâmetros do modelo
# sum(inf$infmat[,"hat"])A partir da tabela acima, temos que:
Observação 9: ponto aberrante;
Observação 10: ponto de alavanca;
Observação 27: ponto influente, aberrante e de alavanca.
Ajuste 3
A seguir, tentaremos melhorar o modelo reduzido anterior, agora aplicando a função log à variável resposta preço e selecionando, novamente, as variáveis explicativas através da técnica AIC. Chamaremos esse ajuste de modelo 3.
Modelo 3
model = lm(log(preco)~imposto+areaC+areaT+imposto:areaT,dados)
summary(model)##
## Call:
## lm(formula = log(preco) ~ imposto + areaC + areaT + imposto:areaT,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.177346 -0.047788 -0.008188 0.057666 0.166648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.942587 0.161863 18.179 9.7e-15 ***
## imposto 0.044811 0.018548 2.416 0.0244 *
## areaC 0.150678 0.072639 2.074 0.0500 *
## areaT 0.002498 0.020138 0.124 0.9024
## imposto:areaT 0.001828 0.002492 0.734 0.4708
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08865 on 22 degrees of freedom
## Multiple R-squared: 0.9211, Adjusted R-squared: 0.9067
## F-statistic: 64.17 on 4 and 22 DF, p-value: 8.266e-12MASS::stepAIC(model)## Start: AIC=-126.37
## log(preco) ~ imposto + areaC + areaT + imposto:areaT
##
## Df Sum of Sq RSS AIC
## - imposto:areaT 1 0.004232 0.17714 -127.72
## <none> 0.17291 -126.37
## - areaC 1 0.033819 0.20673 -123.55
##
## Step: AIC=-127.72
## log(preco) ~ imposto + areaC + areaT
##
## Df Sum of Sq RSS AIC
## <none> 0.17714 -127.72
## - areaT 1 0.016153 0.19330 -127.36
## - areaC 1 0.073064 0.25021 -120.39
## - imposto 1 0.191663 0.36881 -109.92##
## Call:
## lm(formula = log(preco) ~ imposto + areaC + areaT, data = dados)
##
## Coefficients:
## (Intercept) imposto areaC areaT
## 2.83045 0.05562 0.18129 0.01509Pelo critério AIC, nosso modelo 3ficará com as variáveis imposto, área construída e área do terreno.
model3 = lm(log(preco)~imposto+areaC+areaT,dados)
summary(model)##
## Call:
## lm(formula = log(preco) ~ imposto + areaC + areaT + imposto:areaT,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.177346 -0.047788 -0.008188 0.057666 0.166648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.942587 0.161863 18.179 9.7e-15 ***
## imposto 0.044811 0.018548 2.416 0.0244 *
## areaC 0.150678 0.072639 2.074 0.0500 *
## areaT 0.002498 0.020138 0.124 0.9024
## imposto:areaT 0.001828 0.002492 0.734 0.4708
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08865 on 22 degrees of freedom
## Multiple R-squared: 0.9211, Adjusted R-squared: 0.9067
## F-statistic: 64.17 on 4 and 22 DF, p-value: 8.266e-12Diagnóstico 3
Normalidade e Envelope
hist(model3$residuals, probability = TRUE, main = 'Resíduos do Ajuste',
xlab = 'Resíduos', ylab = 'Frequência', col = "gray",)
lines(density(model3$residuals), col = "slategray", lwd = 4)
hnp(model3, main = 'Gráfico Normal com Envelope',
xlab = 'Quantis Teóricos',
ylab = 'Resíduos',
pch = 16)## Gaussian model (lm object)
Análise de Resíduos
plot(lm(model3, dados), which = 1:4, pch = 16,lwd = 3,
caption = c("Resíduos x Valores ajustados",
"Normal Q-Q",
"Resíduos padronizados x Valores ajustados",
"Distância de Cook",
"Resíduos x Alavancagem"))
Medidas de Influência
DFBetas (dfb): estatísticas que indicam o efeito da remoção de cada observação sobre as estimativas dos parâmetros do modelo (apontam os pontos influentes);
DFFit (dffit) e Cook’s D (cook.d): são estatísticas que indicam o efeito da remoção de cada observação sobre os valores preditos/ajustados do modelo (apontam os pontos aberrantes);
COVRATIO (cov.r): estatística que indica o efeito da remoção de cada observação sobre a matriz de covariâncias do modelo, em outras palavras, mede a alteração na precisão das estimativas dos parâmetros do modelo (aponta os pontos alavanca);
HAT (hat): diagonal da matriz de projeção (\(H = X(X'X)^{-1}X'\)) da solução dos mínimos quadrados.
inf = influence.measures(model3)
summary(inf)## Potentially influential observations of
## lm(formula = log(preco) ~ imposto + areaC + areaT, data = dados) :
##
## dfb.1_ dfb.imps dfb.areC dfb.areT dffit cov.r cook.d hat
## 8 0.13 1.04_* -0.04 -1.51_* -1.78_* 0.42_* 0.59 0.26
## 9 0.17 -0.03 -0.17 0.10 -0.28 2.70_* 0.02 0.56_*
## 10 -0.39 0.04 0.10 0.17 0.49 1.74_* 0.06 0.36
## 27 -0.33 -2.17_* 1.56_* 0.68 -2.26_* 2.35_* 1.20_* 0.67_*# A soma da diagonal da matriz H é o número de parâmetros do modelo
# sum(inf$infmat[,"hat"])A partir da tabela acima, temos que:
Observação 8; ponto influente, aberrante e de alavanca;
Observação 9: ponto de alavanca;
Observação 10: ponto de alavanca;
Observação 27: ponto influente, aberrante e de alavanca.
Interpretação
model3$coefficients## (Intercept) imposto areaC areaT
## 2.83045293 0.05562402 0.18129346 0.01509407\[ preço = 2.83 + 0.05*imposto + 0.18*areaC + 0.015*areaT \]
Por fim, com o modelo 3, temos que o preço de venda sobe, em média, 0.05 a cada unidade aumentada na variável imposto, 0.18 na variável área construída e 0.015 na variável área do terreno.